Sedimentary Signatures of Persistent Subglacial Meltwater Drainage From Thwaites Glacier, Antarctica

Introduction

Thwaites Glacier drains directly into the Amundsen Sea embayment (ASE) of West Antarctica and is undergoing rapid changes including grounding-zone retreat, dynamic ice thinning, and increased ice discharge into the ocean (Milillo et al., 2019; Rignot et al., 2019). The glacier catchment extends deep into the interior of the West Antarctic Ice Sheet (WAIS; Holt et al., 2006), indicating destabilization of Thwaites Glacier has the potential to trigger widespread collapse of the WAIS (Weertman, 1974; Schoof, 2007). Increased community efforts are now focused on improving predictability of how Thwaites Glacier will likely respond to continued climatic warming, in particular over its ice shelf and grounding zone (Scambos et al., 2017). While numerous studies focus on ice-shelf and grounding-zone sensitivity to the incursion of warm, Circumpolar Deep Water (e.g., Joughin et al., 2014; Seroussi et al., 2017; Hoffman et al., 2020), less attention has been paid to understanding ice-shelf and grounding-zone behavior in response to subglacial meltwater drainage beneath and at the grounding zone of Thwaites Glacier.

Subglacial meltwater expulsion into the ocean and subsequent buoyancy-driven turbulence is linked to both enhanced ice melt at or near the grounding zone (Gourmelen et al., 2017; Wei et al., 2020; Nakayama et al., 2021) and ice-shelf basal channel formation (Le Brocq et al., 2013; Marsh et al., 2016). The seaward extension of subglacial channels carved upward into the ice-shelf base encourages concentrated melt by warm marine waters (e.g., Alley et al., 2016) and may lead to ice-shelf fracture (Vaughan et al., 2012; Alley et al., 2016; Dow et al., 2018) and ultimately collapse (Borstad et al., 2012; Borstad et al., 2017). Persistent, channelized subglacial drainage at paleo-grounding zones is hypothesized to inhibit sediment deposition needed to construct stabilizing ice-marginal landforms (i.e., moraines and/or grounding zone wedges; Simkins et al., 2017, 2018). Additionally, records from the Ross Sea seafloor suggest grounding zones retreated greater distances when subglacial channels were actively draining, compared to retreat that occurred when channels were inactive (Simkins et al., 2021).

Today, geophysical and remote sensing observations support an active subglacial plumbing network beneath Thwaites Glacier (Schroeder et al., 2013; Smith B.E. et al., 2017; Hoffman et al., 2020). Basal shear stress (Joughin et al., 2009) and regionally elevated geothermal heat flux (Damiani et al., 2014; Dziadek et al., 2021) contribute to meltwater production. Numerically estimated steady-state subglacial water fluxes within the Thwaites catchment range from ∼1 to ∼60 m³/s (Le Brocq et al., 2013). These values are similar to estimates of steady-state subglacial meltwater flux along the Siple Coast that reach the grounding zone through both distributed and channelized drainage pathways (Carter and Fricker, 2012). While no observational evidence of point-source meltwater delivery to the grounding line currently exists, subglacial hydropotential models predict persistent subglacial channels that drain at a few discrete locations along the grounding line (Le Brocq et al., 2013; Hager et al., 2021). Seaward of modern Thwaites Glacier, bathymetric channels connecting paleo-subglacial lake basins indicate subglacial water also flowed beneath the formerly-expanded Thwaites Glacier and neighboring Pine Island Glacier (Lowe and Anderson, 2003; Nitsche et al., 2013; Witus et al., 2014; Kuhn et al., 2017; Kirkham et al., 2019; Hogan et al., 2020). However, the morphology of these bedrock channels is hypothesized to reflect high-magnitude drainage events over several glacial cycles (Kirkham et al., 2019), and therefore do not provide insight into frequency or magnitude of recent drainage events on decadal to millennial timescales.

Subglacially-sourced meltwater plumes mobilize and distribute glaciogenic sediments into the ocean. Sediment cores from the continental shelf can thus be used to provide details of meltwater drainage events that cannot be inferred from observational records or bathymetric landform analysis. Meltwater plume deposits on the continental shelf are characterized as silt-rich sediments with high water content and low biogenic and ice-rafted debris (IRD) content, low shear strength, and low magnetic susceptibility (Kirshner et al., 2012; Hillenbrand et al., 2013; Witus et al., 2014; Smith, J.A. et al., 2017; Streuff et al., 2017; Simkins et al., 2017; Prothro et al., 2018, 2020; Majewski et al., 2020). These properties allow for the distinction of meltwater plume deposits from subglacial, ice-proximal, and seasonally open-marine sediment facies. Incorporation of meltwater plume deposits into traditional facies models of glacier retreat has clarified a connection between subglacial hydrologic activity and grounding-zone position (Simkins et al., 2017, 2021; Prothro et al., 2018). This study employs a suite of sedimentological, geochemical, and statistical analyses of sediment cores from seaward of Thwaites Glacier to identify evidence of subglacial meltwater drainage events into the surrounding ocean. We seek not only to characterize styles of subglacial meltwater drainage that may be distributing sediments offshore of other Antarctic glaciers, but evaluate the relative persistence of meltwater drainage from Thwaites Glacier through the pre-satellite era. These advancements in knowledge will help illuminate connections between meltwater drainage and grounding-zone conditions, thereby broadening the framework within which contemporary subglacial hydrologic observations are interpreted.

Materials and Methods

The Thwaites Offshore Research project of the International Thwaites Glacier Collaboration conducted two cruises aboard the RVIB Nathaniel B. Palmer during the austral summers of 2019 and 2020 (NBP19-02 and NBP20-02). High-resolution multibeam bathymetric and acoustic sub-bottom surveys were conducted on both cruises with a hull-mounted Kongsberg EM122 echo sounder at 12 kHz frequency (Figure 1; Hogan et al., 2020) and a hull-mounted Knudsen Chirp 3260 3.5 kHz sub-bottom profiler, providing geomorphological and acoustic stratigraphic context for sediment cores used in this study. Sediment pockets in bathymetric lows and sediment drapes were targeted for coring. Based on core site locations, core recovery, and on-board core descriptions, five sediment cores recovered with a Kasten corer seaward of the Thwaites Glacier Tongue (TGT) and Eastern Ice Shelf (EIS) were selected to examine spatiotemporal evidence of meltwater drainage and associated sedimentation processes (Figure 1; Table 1). The cores were photographed, lithologically described, and sampled on board.

FIGURE 1. Regional map and locations of sediment cores used in this study. Bathymetry of the eastern Amundsen Sea Embayment (A) and in the vicinity of Thwaites Glacier (B). Data in Figure 1B were primarily collected on cruises NBP19-02 (Hogan et al., 2020) and NBP20-02 and are gridded at 50 m. Bathymetry defining much of the eastern margin of the Eastern Ice Shelf was collected during the iStar project cruise JR294. Additional bathymetry and bed topography are from BedMachine (Morlighem et al., 2020), and sub-ice shelf bathymetry is from Jordan et al. (2020). Red and blue arrows denote inflow of warm water and outflow of meltwater-laden fresher water, respectively, from Wåhlin et al. (2021). Grounding-line data are from Rignot et al. (2011) and Milillo et al. (2019) and coastline data are from the Antarctic Digital Database (ADD). Subglacial flux is from the Quantarctica data package and was derived from a subglacial water routing algorithm run at a 1-km resolution, originally published by Le Brocq et al. (2013). Stars approximate where modeled subglacial channels coincide with the Thwaites Glacier grounding zone (Hager et al., 2021).

TABLE 1. Location, water depth, and recovery of the five Nathaniel B. Palmer (NBP) cores investigated in this study. Bathymetric highs H2 and H1 as defined by Hogan et al. (2020). KC = Kasten core; TGT = Thwaites Glacier Tongue; EIS = Eastern Ice Shelf.

Bulk sediment samples for trace-elemental and grain-size analyses were taken at 10-cm depth intervals using plastic sampling tools from ∼2 cm inward from the core barrel to avoid metal contamination. Archive trays from the cores are stored in the Antarctic Core Collection at the Oregon State University (OSU) Marine Geology Repository. A GEOTEK XZ system at the OSU repository was used to measure magnetic susceptibility (MS) with a point sensor at 1-cm resolution downcore. Computed tomography (CT) scans were acquired through the OSU College of Veterinary Medicine and grayscale images were processed using the SedCT package in Matlab (Reilly et al., 2017). Where CT imagery revealed suspected plume deposits and other distinct sediment types between 10-cm sampling intervals, additional 1-cm thick samples were collected and underwent the same suite of analyses described above.

Sediments were wet-sieved for calcareous foraminifera in 5- to 8-cm intervals. Where present, mixed-assemblage calcareous foraminifera, which are very rare in Antarctic shelf sediments, were utilized for radiocarbon dating with the MICADAS system at ETH Zurich (Wacker et al., 2010). Radiocarbon ages in ¹⁴C years were calibrated using the CALIB 8.2 Marine20 curve (Stuiver et al., 2022) and corrected for marine reservoir effects using the age from live (mixed assemblage) calcareous foraminifera in seafloor sediments (1,117 ¹⁴C years) from the area (Table 2). Collectively, these constraints provide an absolute age framework for which to interpret meltwater drainage persistence.

TABLE 2. Raw and calibrated radiocarbon dates. Calibration was performed using Marine20 in CALIB v.8.2. A marine reservoir age of 1,117¹⁴C years from living (rose bengal stained) mixed-assemblages of foraminifera in surface sediments of core NBP19-02 MC16 was used to obtain a ΔR value of 617 (Figure 1B).

Grain-Size and Smear Slide Analyses

Laser diffraction particle analysis was carried out using a Bettersizer S3 Plus that is equipped with a green laser optimized for fine-grained sediments. Representative subsamples (∼0.5 g wet sediment) were treated with sodium hexametaphosphate and allowed to deflocculate in deionized water for a minimum of 48 h prior to analysis. A magnetic stir bar was used to create a homogenized slurry from which aliquots were pipetted into the Bettersizer S3 Plus reservoir. Each aliquot run consisted of seven individual readings that were averaged and binned by cumulative volume percentages. The maximum (D100) and median (D50) grain size, along with other distribution percentiles, for each sample were identified by Bettersizer S3 Plus software and used to calculate grain-size statistics. Grain-size classifications and statistics follow Folk and Ward (1957). We note that, due to pipette width and instrument limitations, the grain-size measurements do not include material coarser than 3 mm and thus larger grains (i.e., clasts) are not represented in the grain-size distributions. This caveat reduces the efficacy of this method for capturing potential IRD. Interpretations of grain-size distributions are therefore supplemented with CT images, which fundamentally discern material based on density and allow us to comment on relative frequency of downcore IRD (c.f., Cederstrøm et al., 2021).

Prior characterizations of meltwater deposits from Pine Island Bay (Witus et al., 2014) and Ross Sea (Simkins et al., 2017; Prothro et al., 2018) used a Malvern Mastersizer for grain-size analysis. In order to compare those deposits to sediments collected offshore from Thwaites Glacier, duplicate samples from the Ross Sea were run on the Bettersizer S3 Plus (n = 3). Representative smear slide analysis of Thwaites Glacier sediment cores was conducted under plane-polarized light to characterize the general composition of the measured grain-size modes.

Physical and Chemical Analyses

Shear strength was measured with a hand-held shear vane at ca. 10-cm depth intervals on board. Water content was determined on bulk sediment samples (minimum of 3.5 g wet sediment) by oven-drying for at least 48 h at 60°C and calculating the difference in weight percent. Dried samples were then powdered by gently crushing with a ceramic mortar and pestle and sieved through a 125-micron stainless steel mesh. Homogenized samples finer than 125 µm were analyzed for trace and rare Earth elements with x-ray fluorescence (XRF) using a ThermoFisher Scientific Niton XL5 Handheld Analyzer. While hand-held XRF analyzers are unable to measure elemental concentrations on discrete sediment samples with the precision of conventional XRF or other spectrometry methods (Dunnington et al., 2019; Hahn et al., 2019), the results allow comparison of general concentrations downcore and between cores. Furthermore, discrete XRF measurements allow for direct comparison of various measurements conducted on the same sediment samples and are less prone to porewater bias when compared to XRF scanner data. Instrument calibration and performance was tested using a NIST 2709 soil standard prior. The elemental concentrations in parts per million (ppm) were measured by taking the average of three energy spectra filters with a 90-s sampling window for each filter. For those elements used in analysis, two sigma values ranged from ± 8.08% to ± 0.18% (mean ± 2.34%), with the exception of Ni which has two sigma values ranging from ± 27.1% to ± 15.5% in cores KC33 and KC08, respectively (mean ± 18.5%). Discrete MS was measured with a Bartington MS2 dual-frequency sensor on the homogenized samples used for the XRF analyses. Readings were corrected for drift and sample mass, and the final values were used in statistical analysis.

Statistical analysis

To compare sedimentological and geochemical analyses between and within cores through a non-qualitative approach, principal component analysis (PCA) was conducted on all samples using the “prcomp” function in R. PCA is a statistical method useful in identifying trends and the factors driving variability within large, multivariate datasets. Input variables included physical sediment properties (water content, discrete MS), grain-size parameters, and trace element concentrations. All measurements input into this analysis were collected on the same sample or sub-sample at the University of Virginia. Variables were centered and scaled (i.e., normalized) within the code to prevent bias based on variable magnitude, and the package “factoextra” was used to create a biplot visualizing statistical results. Collectively, the results of PCA were interpreted with geochemical and physical sedimentological results to establish sediment lithofacies.

Results

Grain-Size Distributions

Grain-size distributions for all samples reveal striking consistency, with a prominent mode at approximately 5 microns (hereafter referred to as Mode 1) present in nearly all, and dominant in most, samples (Figure 2A). A secondary mode at 11 microns (Mode 2) also occurs in varying concentrations. Mode 2 typically increases at the expense of Mode 1 and comprises the highest volume percent in cores closest to the Thwaites Glacier ice-shelf margin (KC04, KC23; Figure 1B). Additionally, many sediment intervals feature a minor (less than 2% by volume) sub-micron component. Coarse “tails” in the distributions, often interpreted in glaciomarine settings as IRD (e.g., Witus et al., 2014) are observed most commonly in cores closest to the ice-shelf margin and, in few discrete instances, in one relatively distal core (KC08). Overall, cores are dominated by fine silt and clay, and the relative contribution of the sand fraction is minimal and varies with proximity to the modern calving line. The instrument comparison conducted on meltwater deposits from the Ross Sea shows the Bettersizer was able to resolve the previously described 10-micron mode into a primary grain-size mode and secondary grain-size “shoulder” identical in size and volume percent to those observed in this study (Figure 2B). Smear slides show that, in all cores, grain size modes one and two are composed of terrigenous grains, and in-tact or fragmented microfossils are rare or absent except in some surface sediment samples (Figure 3). Observed mineral assemblages vary slightly between cores, but are in all cases dominated by low relief, transparent grains presumed to be quartz (Figure 3).

FIGURE 2. Grain-size distributions of sediment core samples (n = 196) on a semi-log axis plotted against volume percent. Core locations are denoted by color in the inset map in panel (A). (A) All samples exhibit varying degrees of ∼5 micron Mode 1 and ∼11 micron Mode 2, similar to meltwater deposits described elsewhere along the continental shelf in the Pacific Sector of Antarctica. (B) Instrumental comparison for grain-size distributions of previously published meltwater deposits from the western Ross Sea (see inset in Figure 1A) collected during cruises NBP98-01 (KC54 26 cm) and NBP1502 (KC17 170 cm; KC24 79 cm) and published in Prothro et al. (2018) and Simkins et al. (2017). In the lower panel are samples KC33 200–202 cm and KC67 50–52 cm from Thwaites Glacier (TG) shown for reference.

FIGURE 3. Smear slide images and corresponding grain-size distributions for select samples. Images were taken under plane-polarized light. Grain-size distributions shown for corresponding intervals have the same axis labels as in Figure 2B. Panels examine grain-size Modes one and two in intervals (A) within the same core, (B) from cores from different geographic locations, and (C) in surface sediments. (A) Core KC04 is dominated by poorly-sorted, ice-shelf proximal sediments (left) but contains discrete layers of purely terrigenous and relatively well-sorted fine silt (right). (B) KC08 (left; offshore from TGT) and KC23 (right; EIS corner). Grain size modes one and two are caused by terrigenous grains that are present in variable amounts despite geographical difference between sites. (C) KC67 is dominated by Mode-1 grains and features sparse, complete diatom frustules.

Sediment Geochemistry

Downcore geochemistry reflects distinct geographic regions offshore of the eastern and western margins of Thwaites Glacier and variable proximity to the ice-shelf margin (Figure 4). Cores near the EIS are characterized by lower Ti and Sr, and elevated Th and Rb compared to those collected near the TGT. Zr and Rb appear strongly controlled by grain size, as more distal cores with very low sand fractions have higher Rb and lower Zr concentrations while ice-shelf proximal cores have the highest concentrations of Zr. Additionally, more ice-distal deposits have lower Ca and higher Fe concentrations than the ice-shelf proximal cores from the same geographic region. All cores have comparable concentrations of Cr, Ni, and V, trace elements common in intermediate and mafic minerals. KC04 shows a prominent enrichment of Ni at 150 cm that is not mirrored by other geochemical or sedimentological data. Ti content (ppm) in bulk sediment is also plotted against Th/Zr and Sr/Th, ratios of immobile trace elements that are unique for particular geologic settings and not significantly influenced by physical weathering processes (e.g., Hawkesworth et al., 1997; Figures 4B,C). Geographic clusters emerge, where data from cores collected from offshore either the TGT or the EIS plot close together (Figure 4C), as do samples within individual cores (Figures 4B,C). Core KC33 is the only exception, because the geochemical composition of some samples more closely resembles sediments from the EIS cores whereas others are geochemically more similar to that of the TGT cores (Figures 4B,C). While Th/Zr ratios for core KC23 are lower than ice-distal cores from the same region (Figure 4B), this is likely explained by the relatively elevated sand fraction throughout this core (Figure 2A) and a subsequent Zr enhancement (cf., Croudace and Rothwell, 2015).

FIGURE 4. Sediment core geochemistry. (A) Downcore concentrations of selected trace elements in ppm. Core legend is given in the upper-right corner of the Ni plot. (B, C) Immobile trace elemental ratios Th/Zr and Sr/Th plotted against Ti. Element ratios of surface sediment samples (<63 microns) representative for the detrital input into the eastern Amundsen Sea Embayment (ASE, sites 30–33 of Simões Pereira et al., 2018) and the western ASE (sites 34, 36, and 38 of Simões Pereira et al., 2018) in the vicinity of Thwaites Glacier are shown for reference. The vertical dashed line indicates a division between eastern and western geochemical signals.

Downcore Physical Properties and Sedimentological Structures

Core KC08 was collected from a deep trough west of the TGT and is largely composed of rhythmic layers of silt-rich sediments with high water content (max. 49%, mean 40%; Figure 5). Lithologic descriptions suggest the layers are normally graded, which is supported by grain-size distributions that show finer and coarser modes varying in volume percents (Figure 2A). Additionally, this core shows the greatest standard deviation of mean grain size (10.4 µm ± 7.11; Figure 5), reflective of grain-size variation within the graded layers. CT imagery clearly shows these graded silt layers and reveal layer thicknesses range from <0.5 to 8 cm (Figures 6A,B). We do not observe these graded silt layers in other sediment cores analyzed for this study. While overall sand content is low or zero (average 2.9%), two exceptions that elevate mean grain size coincide with the bases of graded silt layers at 25–26 cm and 273–274 cm (11% and 25% sand, respectively; Figure 5). MS increases slightly downcore with occasional excursions and one prominent maximum at 233 cm, coincident with reduced water content within a discrete unit of stratified diamicton (Figure 5). A single radiocarbon date of 520 (2σ 318–674) cal BP (Table 2) comes from an interval with several graded silt layers (52–60 cm; Figures 6A,B). Shear strength increases slightly downcore but does not exceed 2 kPa.

FIGURE 5. Downcore lithology, physical sediment properties, and sedimentary structures for each core. GS = grain size; MS = magnetic susceptibility; WC = water content; SS = shear strength. Sedimentary structures and stratigraphy are further discussed in Figure 6 with the added context from computed tomography images. Linescan MS profiles and discrete MS profiles are plotted together as a line and closed circles, respectively. Discontinuous linescan MS profiles can result from uneven core surfaces (e.g., KC04). Note the difference in scale of MS and Mean GS between cores.

FIGURE 6. Greyscale computed tomography (CT) images of sediment cores collected offshore Thwaites Glacier. (A) CT image next to linescan photo. Downcore Sr/Th is imposed on core KC33 and reflects a gradual change in sediment source. Dates (cal yr BP) are shown alongside the corresponding interval. Refer to Table 2 for error bounds. Colored vertical lines denote core sections shown in (B–F). (B–F) Sediment facies described in this study with annotations. Long-axis orientation for clasts of varying sizes is traced in diamicton facies (C, F). Note that more than one sediment facies is observed within each core.

Core KC04, collected on bathymetric high H2 ∼20 km from KC08, exhibits high variability in most measured properties (Figure 5). Water content is ∼41% in the upper 30 cm and decreases to ∼26% downcore with the exception of two discrete intervals (88–89 cm, 210–212 cm) where water content increases to ∼35%. MS fluctuates downcore, reflective of the highly variable sand fraction (range: <1–34%), while mean grain size shows little downcore variability (13.7 µm ± 4.38). This matrix material is poorly sorted (mean 1.84) with a maximum of 3.21 at 46–47 cm (very poorly sorted). Core stratigraphy reveals several sediment types including massive and stratified diamictons (Figure 6A). Fine-grained, parallel to sub-parallel laminations and beds ranging from <1 to ∼10 cm in thickness are observed within stratified diamicton, as are gravel and pebbles oriented with their long axes parallel to lamination (Figure 6C). Pebbles largely disappear upcore from the most poorly-sorted interval at 46–47 cm. Dates from three intervals between 130 and 62 cm core depth span 490 (2σ 302–654) to 270 (2σ 26–478) cal BP and are in stratigraphic order (Table 2). Shear strength increases slightly downcore from ∼0.25 kPa to 2 kPa. Core KC33, collected farthest offshore of the EIS, is characterized by high water content (range, 37–52%; mean 46%) and low shear strength that increases slightly downcore, coincident with a modest decrease in water content (Figure 5). The mean grain size hardly varies downcore (7.37 µm ± 1.77) although volume percent of Mode 2 grains increases slightly upcore from ∼50 cm (Figures 2A, 5). Sediments are, on average, better sorted than the matrix material of ice-proximal diamictons (mean 1.29). The core base shows faint, irregular lamination that transitions into sediment which appears bioturbated and mottled from ∼210 cm core depth to the core top (Figure 6D). MS varies moderately, but remains generally low except for peaks at 217 and 257 cm that align with coarse-grained layers apparent in the CT scan (Figure 6A). Calcareous microfossils from the core base (290–295 cm) yield a single date of 9,630 cal BP (2σ 9,325–9,967).

Sediments in core KC67 from the eastern side of the EIS have high water content (range, 42–51%; mean, 46%), similar to cores KC33 and KC08, and consistently low MS (Figure 5). Shear strength fluctuates between 0 and 2.5 kPa and increases slightly downcore, while the average grain size remains remarkably uniform throughout (5.53 µm ± 0.28). Sediments show strong laminations with sharp contacts, interspersed with scarce IRD below ∼130 cm core depth (Figure 6E). Above this depth, sediments are heavily bioturbated. The clear lamination suggests alternations between thin coarser- and finer-grained layers, but the 1–2-cm thickness (i.e., depth range) of the grain-size samples is too large to resolve grain-size differences between individual laminae (Figure 2A). Overall, KC67 exhibits the lowest (i.e., most well-sorted) sorting values that are consistent downcore (mean 1.09; Figure 5).

East of the EIS, core KC23 consists predominantly of massive diamicton (Figures 5, 6A,F). The diamicton below 15 cm core depth is characterized by low water content (mean 18%), fluctuating MS, and a strikingly consistent grain-size matrix (11.5 µm ± 1.29; Figure 5). Mean sorting of the diamicton (1.81) is very similar to the matrix material of KC04. Between ∼8 and 15 cm core depth, a normally-graded, oblique sand bed ∼5 cm thick separates the diamicton from overlying silty-clays with notably higher water content (average 39%). MS is lower in the silty-clay unit than in the sand bed and the diamicton. Shear strength through the diamicton is low (≤1 kPa) but increases up to 4 kPa below 140 cm (Figure 5).

Principal Component Analysis

The first two principal components (PC1 and PC2) account for 40.4% and 27.1% of variability in the dataset, respectively. Samples are plotted along PC1 and PC2 axes, along which sample positioning is controlled by the magnitude and relative sign (i.e., positive or negative) of the eigenvectors driving the variance (Figure 7). For PC1, the eigenvectors for Fe, V, water content, and clay percent push samples into negative PCA space, and sand percent, D100, D50, and Ca drive samples into positive PCA space. Negative space of PC2 is largely influenced by Rb, Th, Ti, and clay and silt percentages, while positive PC2 space is defined by Sr, MS, and Zr. Grain-size parameters appear to be the dominant control over PC1, separating finer-grained samples with higher water content, total Fe, and clay-sized grains from samples with higher sand fraction, more variable sorting, and overall coarser grains. While grain-size parameters also influence PC2, geochemistry seems to play a more important role as cores are divided by geographic region and presumably by sediment source. Cores KC08 and KC04 straddle the zero line of PC1, indicating these cores contain intervals deposited by variable sedimentary processes. In contrast, KC33 and KC67 plot entirely in negative PC1 space with the exception of two sample outliers in KC33. Core KC23 samples cluster tightly except for sample outliers including the surface sediments (Figure 7, denoted by an enlarged icon) that plot in negative PC1 space and others that plot high on the PC1 axis. These outliers correspond to samples collected from the sorted sand bed (Figures 5, 6A).

FIGURE 7. Biplot projecting results of PCA. Variance in principal component (PC) one is largely controlled by grain-size parameters and water content while PC2 variance is influenced by trace element concentrations. The surface sample (0–2 cm) for each core is indicated by an enlarged symbol, and concentration ellipses are drawn to include the 95th percentile.

Facies Interpretation

The multi-proxy analysis employed in this study, including CT imagery and the integration of PCA, is valuable in characterizing and distinguishing between sediment types recovered offshore of Thwaites Glacier. We interpret PC1 as distinguishing glacially-sourced sediment (i.e., via ice rafting, basal melt out, gravity flows, or till from the glacier base) from sediment transported by subglacial meltwater processes, including within plumes emanating from the grounding line. High water contents, clay percentages, and Fe contents are shown to be statistically important characteristics for identifying meltwater plume deposits (Figure 7), largely consistent with studies characterizing such deposits from the Antarctic continental shelf (e.g., Witus et al., 2014; Simkins et al., 2017; Prothro et al., 2018). Sediment samples from KC33 and KC67 plot almost exclusively in negative PC1 space (Figure 7), and thus we interpret these cores to reflect pure meltwater endmembers that manifest as both laminated and bioturbated sediment facies.

Two diamicton facies, massive and stratified, were recovered in ice-shelf proximal cores KC04 and KC23. Additionally, a discrete unit of stratified diamicton is observed in the lower part of KC08 (Figure 6A). Clast orientation in the stratified diamicton, along with observed soft-sediment deformation and sediment draping above clasts, is characteristic of ice rafting, where debris released from the base of floating ice settles through the water column. The sediments comprising laminations and beds in stratified diamicton are rich in Mode-1 grains (Figure 3A). In glaciomarine settings, diamictons indicate deposition subglacially or through processes at or near the grounding zone (e.g., Domack and Harris, 1998; Prothro et al., 2018; Smith et al., 2019). Because this study is focused on sediments that inform meltwater processes, we do not explore the specific depositional processes responsible for recovered diamictons in depth.

Much of the sediment from KC08 plots in negative PC1 space (Figure 7) and shares similarities with sediments from KC33 and KC67, including high water contents and high Fe contents (Figure 5). We therefore interpret core KC08 as consisting of primarily meltwater-derived sediments, composing a normally graded silt facies that reflects deposition from hyperpycnal flows, possibly as turbidity currents. To further characterize this unique facies offshore of the TGT, we examine the downcore frequency of different sediment types found in cores KC08 and KC04. Within units of stratified diamicton, both beds and individual, continuous laminae are counted collectively as “laminations”because of the similarity in grain-size. Laminae and normally-graded silt layers are counted in 10-cm intervals, and layers that cross an interval boundary are counted in both intervals (Figure 8A). The abundance of graded silt layers within KC08 fluctuates slightly but remains common throughout (mean 4.7 layers per 10-cm depth intervals), while few instances of similar layers are observed in the upper 120 cm of KC04 (Figure 8A). We also quantify downcore variability in thickness of graded silt layers in KC08, and find that 75% of the layers are 4 cm thick or less, while outliers with thicknesses between 4.5 and 8 cm thick are concentrated in the upper 70 cm of the core (Figure 8C).

FIGURE 8. Frequency analysis of sediment layers observed in cores offshore the TGT. (A) Downcore frequency of sediment facies categorized as massive diamicton (Dm), graded silt layers and stratified diamicton (Ds), in addition to laminations and beds within facies Ds (lams). (B) Representative grain-size distributions for each sediment facies. Note that ice-rafted debris in Ds includes only the sand fraction and not coarser clasts. (C) Thickness of graded layers throughout core KC08. Boxplot denotes 75% of all graded layers are ≥2.5 cm. While thick (i.e., ≥5 cm) graded layers are concentrated in the upper 70 cm of the core, a comparable layer is observed at ∼280 cm depth.

Discussion

Surface sediments throughout Pine Island Bay, which receives drainage from both Thwaites and Pine Island glaciers, are dominated by meltwater plume deposits (Witus et al., 2014), indicating the sustained and widespread importance of sediment delivered by meltwater drainage into the ASE. Meltwater deposits are particularly thick in basins and deep channels on the inner continental shelf (Lowe and Anderson, 2003; Nitsche et al., 2013; Kuhn et al., 2017), and become thicker with increasing proximity to the modern grounding zone of Pine Island Glacier (Witus et al., 2014). This spatial pattern of meltwater deposits indicates uniform dispersal by sediment-laden plumes, and accumulation of plume deposits in bathymetric lows, possibly remobilized across the seafloor by hyperpycnal flows or downslope gravity flow processes (e.g., turbidity currents). In Marguerite Bay (Kennedy and Anderson, 1989) and the in Ross Sea (Simkins et al., 2017; Prothro et al., 2018), similar meltwater deposits overlie proximal glaciomarine sediments and are capped by or interspersed with distal, diatom-bearing glaciomarine muds. This stratigraphic sequence is characteristic of the most recent deglaciation of the Antarctic continental shelf and onset of open-marine conditions, typically as the result of rapid grounding-zone retreat (Bentley et al., 2011; Bart et al., 2017; Prothro et al., 2018). However, unique to Pine Island Bay is the absence of the diatomaceous mud which is not observed in any cores described to date (Lowe and Anderson, 2003; Kirshner et al., 2012; Hillenbrand et al., 2013; Nitsche et al., 2013; Smith et al., 2014; Witus et al., 2014; Kuhn et al., 2017), and is similarly absent in the cores used in this study. Instead, surface sediments from each core analyzed here are characterized by high water content (Figure 5), a prominent Mode 1 (Figure 2A), and plot in negative PC1 space (Figure 7), indicating that meltwater drainage is the primary mode of sedimentation seaward of contemporary Thwaites Glacier.

Here, we discuss the source of meltwater plume deposits offshore Thwaites Glacier, the variability in styles of meltwater drainage, and consider processes associated with meltwater evacuation, as well as floating ice and oceanographic responses. Calibrated radiocarbon dates presented here (Table 2) are useful in conjunction with core stratigraphy for evaluating persistence and frequency of subglacial meltwater drainage through time (Figure 6A). Radiocarbon at the base of core KC33, the northern-most core in this study, suggests sediment accumulation on the H1 high spanning the last ∼9,000 to 10,000 years (Table 2). This is consistent with regional reconstructions of deglaciation (Larter et al., 2014) and radiocarbon dates that indicate grounding-zone retreat from the H1 high by ∼10,300 years ago (Hillenbrand et al., 2017; Hogan et al., 2020). In the vicinity of Thwaites Glacier, sedimentation rates for similar glaciomarine facies range from 0.095 cm/a beneath the Pine Island Ice Shelf (Smith et al., 2019 and references therein) to 0.02–0.086 cm/a for core NBP99-02 TC49 which sampled meltwater plume deposits in inner Pine Island Bay (Figure 1B; Witus et al., 2014). Altogether, these data allow for meaningful considerations of sediment depositional processes through time.

Source and Grain-Size Production of Meltwater Deposits

Rates of sub-ice shelf melt from observations (17.7 m/yr; Rignot et al., 2013) and models (11.6 m/yr; Pelle et al., 2019) of Thwaites Glacier implicate the formation of buoyant, meltwater plumes without a subglacial origin. Such plumes would, similarly, influence stratification and ocean circulation in an ice-shelf cavity and could entrain and transport debris melted out from the ice-shelf base. Albeit dependent on basal debris thickness and ice-shelf melt rates, debris advected by basal ice is largely melted out within the first few kilometers of the grounding line (e.g., Prothro et al., 2018), restricting the area over which sediment entrainment by ice-shelf melt plumes would be expected. Downstream refreezing would decrease the carrying capacity of sediment-laden plumes and promote deposition beneath the ice shelf. Furthermore, the rate at which a plume entrains ambient water in an ice-shelf cavity (in this instance, warm CDW at depth) is a function of plume speed (Jenkins, 2011). Discharge of subglacial meltwater, therefore, has an inherently greater capacity to drive sub-ice shelf melt through entrainment of warm water than plumes that originate in ice-shelf cavities. Given the widespread distribution of meltwater plume deposits and the distance of meltwater endmember core sites from the recent grounding-line and ice shelf-margin positions (Figures 1B; Supplementary Figure S2), a subglacial source better explains these silt-rich deposits than a basal-ice shelf source. We acknowledge, however, that Mode-1 grains mobilized in sub-ice shelf plumes and deposited some distance offshore Thwaites Glacier would have the geochemical signatures as plume deposits sourced from the subglacial environment.

The unique geochemical signals observed in sediment cores offshore of the EIS and TGT reflect spatial variations in source material, and therefore subglacial geology, beneath the eastern and western portions of Thwaites Glacier (Figure 4). This is consistent with geochemical fingerprinting of fine-grained surface sediments from the inner shelf of Pine Island Bay (Simões Pereira et al., 2018). Planar, seemingly non-erosional contacts and the absence of sand-rich, winnowed deposits in ice-distal cores suggest marine currents have not remobilized sediments deposited at these core sites (Figure 6). Similarly, geochemical consistencies between diamictons recovered near the ice-shelf margin and meltwater plume deposits further offshore (Figure 4) imply Mode-1 rich sediments offshore Thwaites Glacier were not delivered to core sites by currents originating offshore another glacial system. The ubiquity of grain-size Modes 1 and 2, despite spatial variations in subglacial geology, strongly supports that grain-size production of meltwater deposits is a process-controlled, rather than a source-controlled, phenomenon. Previous studies have inferred subglacial till as the source of meltwater plume deposits based on a shared grain-size mode (e.g., Witus et al., 2014; Simkins et al., 2017; Prothro et al., 2018), a relationship that is similarly evident from grain-size distributions of ice-shelf proximal diamicton and ice-distal sediments in this study (Figure 2A). Grain surface micromorphology of meltwater deposits from Pine Island Bay shows evidence of glacial influence without significant surface alteration by hydrologic transport (Witus et al., 2014), suggesting the grain-size distribution of meltwater deposits reflects the process of glacial abrasion rather than size reduction through meltwater transport. The sedimentological likeness of Thwaites Glacier meltwater deposits to those produced by other Antarctic glaciers and paleo-ice streams across a range of space and time implies widespread and consistent subglacial erosion and hydrologic transport to the Thwaites Glacier grounding zone.

Spatial Variability in Subglacial Meltwater Plume Deposition

The stratigraphy revealed by CT imagery suggests contrasting depositional processes associated with subglacial meltwater drainage offshore from the TGT and EIS (Figure 9). Cyclical, normally-graded silt layers dominated by meltwater deposits indicate persistence of sedimentary processes largely unique to core site KC08 (Figure 8A). Calcareous benthic foraminifera are found in several discrete intervals of KC08, but rarely in sufficient numbers for radiocarbon dating. The corresponding layers could reflect infrequent gravitational downslope sediment transport (e.g., Smith et al., 2009) as the core site lies approximately 400 m below the local carbonate compensation depth (Hillenbrand et al., 2003; Majewski, 2013). Alternatively, preservation of calcareous foraminifera could result from rapid sediment deposition that limited exposure to corrosive bottom waters. Contacts separating normally-graded silt layers are sharp but appear non-erosive, suggesting insufficient depositional energy to erode underlying sediments. Additionally, intervals with continuous bands of IRD (Figure 6A) possibly mark periods of meltwater drainage cessation when sedimentation was temporarily dominated by IRD deposition (e.g., Prothro et al., 2018). These combined observations support the interpretation that the normally-graded silt layers in core KC08 were deposited by hyperpycnal flows generated during episodic meltwater drainage events of varying magnitudes (Figure 9B). If we hypothetically consider these silt layers in the context of a fixed grounding-zone position, which can be inferred from the presence of a grounding zone wedge on the northern crest of H2 (Hogan et al., 2020), layer thickness between 0.5 and 4 cm may reflect punctuated, lower-magnitude meltwater discharge, as this describes 75% of such units (Figure 8C). Hence, the thicker (4.5–8 cm) normally graded silt layers may result from relatively higher-magnitude drainage events. In the context of a grounding zone that has been retreating landward, thicker graded silt layers still implicate higher-magnitude drainage events, as in this scenario the source of drainage is moving progressively further from the core site.

FIGURE 9. Cartoon of contrasting meltwater drainage styles inferred offshore from the EIS and TGT. Callouts show general depositional settings, not exact locations, of sediment facies presented in Figure 6. Grounding-line positions are not exact. (A) Persistent, low-magnitude drainage emanating from beneath eastern Thwaites Glacier. Settling from suspension dominates here, likely beneath a perennial canopy of floating ice. (B) Downslope gravitational transport sourced from sediment-laden meltwater plume deposits. Hyperpycnal flow delivers coarse silts and occasionally fine sands to trough, where fine silts and clays settle from suspension.

Sediments in core KC04, recovered from within a ∼60 m deep channel on the H2 high (Figure 1), record the recent unpinning of the TGT (Rignot et al., 2011; Hogan et al., 2020) and the transition to seasonally open-marine conditions at this core site. The stratified diamictons comprising most of the sediments in this core (Figures 6A, 8A) feature laminated and stratified intervals with large clasts and a matrix dominated by Mode 1, and are consistent with deposition in an ice-proximal environment (e.g., Prothro et al., 2018; Smith et al., 2019). Estimated sedimentation rates of stratified diamicton in core KC04 (between 62 and 130 cm core depth; Table 2; Figure 6A) are an order of magnitude higher than sub-ice shelf and grounding-zone distal sediment facies observed in the region (Witus et al., 2014; Smith et al., 2019). Clast orientation is parallel to bedding (Figure 6C), which along with observed soft-sediment deformation and sediment draping above clasts is characteristic of IRD released from the base of floating ice settling through the water column. The prevalence of concentrated meltwater silt laminae within clast-rich diamicton in core KC04 may record subglacial meltwater drainage beneath this region of western Thwaites Glacier when the grounding zone was nearer to the core site (Figure 9B). Lithology transitions upcore of ∼30 cm (Figures 5, 6A) from stratified diamicton to laminated and bedded sediments with IRD-rich layers. These sediments are characterized by higher water content, finer mean grain size, and better sorting compared to the diamicton downcore (Figure 5), all of which are consistent with meltwater plume deposition.

Pervasive bioturbation in meltwater endmember core KC33 is consistent with deposition in an open-marine setting. Various configurations of the TGT over the last ∼12 years, prior to calving of iceberg B22, confirm open-marine conditions at the core site (cf., MacGregor et al., 2012; Supplementary Figure S3). Interestingly, the source of detritus in core KC33 appears to shift upcore from material originating from beneath eastern Thwaites Glacier to material from below western Thwaites Glacier (Figures 4, 6A). This geochemical shift is not mirrored by substantial changes in grain size or other physical properties (Figure 5), suggesting consistent depositional processes for the sediments recovered at this site. Sediment mixing by bioturbation does not explain this shift, as the preserved laminae at the core base and the overlying bioturbated sediments have a distinct geochemical composition (mean Sr/Th from 250 cm core depth to core base = 6.4) from bioturbated sediments at the core top (mean Sr/Th from core top to 52 cm core depth = 8.3; Figure 6A). Given the strong meltwater signature of the sediments, we interpret the shift in geochemical signature to indicate rerouting of subglacial drainage beneath or emanating from Thwaites Glacier. Subglacial channel reorganization may be driven by changes in subglacial hydropotential, resulting from changes in glacier geometry (e.g., Carter et al., 2013) or supply of basal meltwater to the grounding zone (Simkins et al., 2021). Whether this shift in source material in KC33 reflects a migration or shutdown of a former meltwater drainage pathway, or the development of a new preferential drainage pathway cannot be constrained. A change in sediment source could also result from blocking or siphoning of sediment-laden meltwater plumes by the rugged bathymetry on the inner shelf seaward of Thwaites Glacier. Mean grain size increases slightly in the upper ∼40 cm of the core, which may suggest strengthening of ocean currents, but this grain-size change doesn’t occur contemporarily with the geochemical change. Stronger currents transporting meltwater plume deposits from offshore western Thwaites Glacier to the H1 high therefore does not explain this observed change in source material. These findings from KC33 provide evidence for a dynamic subglacial hydrologic network beneath Thwaites Glacier over the last ∼10,000 years and implicate drainage reorganization in association with major ice-configuration changes.

Structurally, the laminations defining over half of the sediments in core KC67, and the lowermost part of core KC33, resemble those of modern sub-ice shelf sediments recovered from below the Pine Island Glacier ice shelf (Smith, J.A. et al., 2017), paleo-sub ice shelf sediments in the Ross Sea (Simkins et al., 2017; Prothro et al., 2018), and sediments offshore the eastern Antarctic Peninsula (Subt et al., 2017) and of Ryder Glacier in Greenland (O’Regan et al., 2021). The consistency of laminae thicknesses in KC67 (Figures 6A,E) suggests a pulsing outflow of sediment-laden subglacial meltwater while the uniform downcore geochemistry (Figure 4) implies drainage along a persistent subglacial pathway. Hydropotential calculations predict regionally-significant fluxes of subglacial meltwater to the grounding zone ∼50 km from the site of core KC67 (Figure 1B; Le Brocq et al., 2013). This suggests that lateral plume migration and interactions between sediment-laden plumes and ocean currents may exert important controls on meltwater deposit distribution offshore the EIS. Preservation of laminations in KC67 and lowermost KC33, which may be destroyed by bioturbation following ice-shelf retreat (Jung et al., 2019), supports deposition beneath a perennial ice canopy–floating ice that may include an ice shelf, ice mélange, and/or fast ice (c.f., Milliken et al., 2009; Minzoni et al., 2015; Lamping et al., 2020). No age constraints are available for core KC67, although nearby core NBP99-02 TC49, which consists of similar meltwater plume deposits (Figure 1B; Unit one from Witus et al., 2014) that accumulated at ∼0.086 cm/yr. While the grounding-zone position is not assumed to be static through the deposition of meltwater-derived sediments in KC33 and KC67, the absence of diamictons and general lack of IRD indicate a depositional environment consistently seaward of the sub-ice shelf basal debris zone (Prothro et al., 2018; Smith et al., 2019). The transition from the laminated to bioturbated facies above ∼130 cm core depth occurs without major changes in grain size (Figure 5), suggesting suspension settling from buoyant meltwater plumes remains the primary mechanism of sedimentation (Figure 9A). Individual burrows are observed within laminated sediments as deep as 190 cm in this core (Figure 6A) which might indicate the presence of infaunal benthic organisms coeval with meltwater drainage beneath a floating ice canopy.

The inferred spatial variation in subglacial drainage style is intriguing, as floating ice and grounding-zone characteristics of the eastern and western parts of Thwaites Glacier also exhibit unique, albeit connected, patterns of behavior (Schroeder et al., 2016; Milillo et al., 2019; Miles et al., 2020; Alley et al., 2021). Overall, meltwater drainage delivering detritus offshore the slow-flowing EIS is characterized as low-magnitude, yet persistent, with important implications for nutrient supply to marine ecosystems (Gerringa et al., 2012; Vick-Majors et al., 2020). Offshore the TGT, extending from the fast-flowing central axis of the glacier, stratified diamictons and normally-graded silts with occasional occurrence of IRD-bearing layers implicate meltwater drainage preceded floating ice-tongue (or ancestral ice shelf) retreat and a relatively recent increase in intensity of meltwater discharge events and volume of freshwater in the ice-shelf cavity.

Connection of Meltwater Plumes to the Contemporary Ice-Ocean System

Release of subglacial meltwater into the ocean is thought to initiate hyperpycnal flows through rapid sediment loading (e.g., Noormets et al., 2009; Sutherland et al., 2020; King et al., 2022). Although such interpretations are largely based on findings from the continental shelf edge and uppermost continental slope, the relative relief of H2, the established position of grounded ice on this high, and presence of sediment fans on its northern flank (Hogan et al., 2020) suggests similar processes are feasible close to Thwaites Glacier (cf., Dowdeswell et al., 2015). In flume experiments, hyperpycnal currents initiated by convective sedimentation (i.e., driven by suspension settling) are most likely to form when the buoyant inflow is carrying a greater particle load (Davarpanah Jazi and Wells, 2020). Conceptually, the concentration of suspended detritus increases with discharge and hydraulic pressure gradients in subglacial channels (e.g., Alley et al., 1997; Swift et al., 2005; Gimbert et al., 2016), thus supporting a possible connection between observed thicker normally-graded silt layers and relatively higher magnitude of drainage events with heightened particle load mobilized by plumes. The intervals of KC08 that bear few calcareous microfossils, including foraminifera, suggest processes including meltwater-induced hyperpycnal transport, recycling and reworking by glacial advance and subsequent transport as IRD or buoyant interflows (i.e., plumes) are all possible. The single age of 520 cal BP (2σ 318–674) from 52 to 60 cm core depth therefore provides a maximum depositional age, and the concentration of relatively thicker graded silt layers above this depth may imply an increase in magnitude of drainage events in recent centuries.

Based on the spatial distribution of core sites and collective data presented here, we find evidence that plume deposits blanket the upper 5–10 cm of seafloor along the EIS and TGT. Thickness of deposits increases substantially (i.e., 3 + meters; the thickness of meltwater endmember cores) with increasing distance from the ice margin. These results suggest that contemporary plumes transport suspended detritus between 5 and 66 km from grounded ice to these core sites (Supplementary Figure S2). Plume deposits on shallow topographic highs (<400 m water depth) and in deep troughs (>1,000 m water depth) indicate plumes travel at relatively shallow depths in the water column, consistent with observations of the glacial meltwater-enriched layer in the ASE (Biddle et al., 2019). Furthermore, recent oceanographic work has detected fresher, glacial meltwater emanating from beneath the TGT (Wåhlin et al., 2021). This water is routed seaward toward the trough from which KC08 was recovered (Figure 1B; Wåhlin et al., 2021). While the source is presumed to be sub-ice shelf and/or grounding-zone ice melt, the sedimentary evidence of meltwater-derived deposition in cores KC04 and KC08 suggests subglacial drainage is an important contributor to this freshwater signal. Conversely, the substantial inflow of warm water to the east of H2 (Figure 1B; Wåhlin et al., 2021) supports the likelihood of enhanced turbulent mixing with fresh meltwater outflow along the grounding zone of western Thwaites Glacier, as is recently suggested for neighboring Pine Island Glacier (Nakayama et al., 2021). Numerical modelling of channelized drainage pathways beneath Thwaites Glacier predicts persistent drainage that discharges beneath the TGT (Hager et al., 2021). In some instances, these predicted locations of channelized drainage coinciding with grounding-zone embayments, which has also been observed in bathymetric records from the Ross Sea (Simkins et al., 2018).

Release of large quantities of freshwater and sediment into the marine system from paleo-subglacial basins can occur due to changes in subglacial hydraulic conditions (Lowe and Anderson, 2003; Nitsche et al., 2013; Witus et al., 2014; Kuhn et al., 2017; Simkins et al., 2017; Prothro et al., 2018). Bathymetry in the modern Thwaites ice-shelf cavity inferred from airborne gravity measurements reveals isolated basins (Millan et al., 2017; Jordan et al., 2020) similar in scale to those seaward of the grounding zone (Nitsche et al., 2013; Hogan et al., 2020) in which subglacial meltwater is likely to have ponded (Witus et al., 2014; Kuhn et al., 2017; Kirkham et al., 2019). It is possible that these paleo-subglacial basins now beneath the contemporary ice tongue served as reservoirs for subglacial meltwater and sediment, and that higher-magnitude drainage events responsible for thicker graded silt packages were triggered by grounding-zone retreat. This model is also called on to explain the extensive meltwater facies observed seaward of Pine Island Glacier (Witus et al., 2014) and in association with previously expanded ice streams elsewhere, such as in Marguerite Bay (Anderson and Fretwell, 2008). The role of basins as reservoirs of subglacial meltwater helps resolve discrepancies between higher-magnitude drainage events implied by the sediment record and annual estimates of subglacial meltwater production within the Thwaites Glacier catchment (Joughin et al., 2009; Le Brocq et al., 2013; Kirkham et al., 2019).

Subglacial lakes and channelized and distributed drainage beneath Thwaites Glacier have been inferred from airborne and on-ice geophysical investigations (Schroeder et al., 2013, Schroeder et al., 2016; Clyne et al., 2020) and remote sensing data (Smith B.E. et al., 2017; Hoffman et al., 2020). We can envision a subglacial plumbing network wherein efficiency and connectivity increase downstream with increasing basal water and ice-surface slope (cf., Schroeder et al., 2013). With the integration of sediment core interpretations, we hypothesize these “higher order” channels transport detritus characterized by Mode 1 in suspension to the grounding zone and release it into the marine environment within discrete meltwater plumes. Channelized subglacial drainage could also be responsible for the distinct and consistent geochemical signal of meltwater plume deposits observed offshore from the EIS and the TGT (Figure 4). Such a conceptual model is in line, both in scale and in space, with outputs from recent hydrologic modeling studies predicting subglacial meltwater pathways and channelized drainage beneath past and modern Thwaites Glacier (Kirkham et al., 2019; Hager et al., 2021). While the “cascading” of subglacial water between lake basins upstream of the modern grounding zone (Smith B.E. et al., 2017) is not associated with significant changes in ice-flow dynamics (Hoffman et al., 2020), the release of such volumes of water to the grounding zone, where it could induce a dramatically different ice response, has not yet been observed. Catchment-wide estimates of basal melt production beneath Thwaites Glacier in its modern (Joughin et al., 2009) and Last Glacial Maximum (Kirkham et al., 2019) configurations vary between ∼5 and 12 km³/a; subglacial meltwater pooling and episodic evacuation from such reservoirs may be responsible for triggering downslope gravity flows and producing the thicker graded silt layers recovered offshore the TGT. As grounding-zone retreat of Thwaites Glacier is expected to continue through the coming decades (Joughin et al., 2014; Seroussi et al., 2017; Yu et al., 2018), changes in subglacial hydraulic potential may cause the release of increasing amounts of freshwater from upstream reservoirs, leading to potential amplification of existing drivers of submarine ice melt and delivery of sediments to the ocean.

Conclusion

Sedimentary evidence seaward of Thwaites Glacier confirms the existence and persistence of a frequently active and historically long-lived subglacial hydrologic system. We find quantitative and qualitative evidence that sediments transported in, and deposited from, meltwater plumes generated beneath Thwaites Glacier are not only similar to those extinct sedimentary systems described offshore from other Antarctic glacier systems, but are also suggestive of subglacial erosive and/or hydrologic processes that occur independent of subglacial geology. Geochemical signatures demonstrate the utility in evaluating meltwater drainage persistence, or transience, of subglacial drainage pathways through time. We identify sedimentary evidence that suggests an increase in magnitude of meltwater drainage beneath the TGT in recent centuries, and speculate that subglacial drainage may have preceded or accompanied recent ice-tongue unpinning. Drainage patterns inferred in association with Western Thwaites Glacier differ from those of Eastern Thwaites Glacier, where persistent, low-magnitude meltwater drainage is long-lasting, yet transient in its preferred pathways. The evidence of persistent drainage beneath floating ice underscores the potential for turbulent mixing at the grounding zone to be an important factor for the Thwaites Glacier system. Omission of processes related to subglacial meltwater drainage from models of grounding-zone retreat (Yu et al., 2018), calculations of ice-shelf melt rates (Seroussi et al., 2017), evaluations of ice-shelf stability (Alley et al., 2021), and flux of freshwater into the Southern Ocean (Adusumilli et al., 2020) likely leads to inaccurate estimations and results. In the context of several prior studies presenting evidence of active subglacial drainage accompanying and preceding grounding-zone retreat, our results demonstrating persistent drainage of subglacial meltwater in recent millennia are potentially another ominous signal pointing toward the impending retreat of Thwaites Glacier.

Data Availability Statement

The datasets generated for this study, including sediment physical property, grain-size, and geochemical data, are available as spreadsheets which are archived by the United States Antarctic Program Data Center, an online data repository for Antarctic projects funded by the National Science Foundation hosted by the Lamont-Doherty Earth Observatory of Columbia University (doi: 10.15784/601514). Supplemental figures can be downloaded as a separate document.

Author Contributions

APL and LS conceptualized the project and designed the methodology. All authors participated in data collection at sea and/or data curation during post-cruise sampling. APL conducted formal analysis and interpreted results along with LS and JA. LW performed radiocarbon measurements. APL, LS, JA, RC, JW, C-DH, JS, AAL, RT, RL, KH, FN, and AG were involved in discussion of results. APL produced figures and drafted the manuscript with significant contributions from LS. JA, C-DH, JS, JW, RL, and KH provided critical feedback and constructive comments. All authors read and approved of the submitted manuscript.

Funding

This work is from the THOR project, a component of the International Thwaites Glacier Collaboration (ITGC). Support from National Science Foundation (NSF: Grant 1738942) and Natural Environment Research Council (NERC: Grant nos. NE/S006664/1 and NE/S006672/1). Logistics provided by NSF-U.S. Antarctic Program and NERC-British Antarctic Survey. ITGC Contribution No. ITGC-072. Internal grants from the University of Virginia to APL and LS also helped support this research.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank the captains and crews of the RV/IB Nathaniel B. Palmer, and the science parties from cruises NBP19-02 and NBP20-02. Many thanks to V. Stanley and the Oregon State University Marine and Geology Repository for their assistance with additional sampling and collection of CT images. We thank D. X. Buskard for her assistance with grain-size data collection. Figures make use of the perceptually-uniform colormap “batlow” developed by F. Crameri (2018). We thank editor B. J. Davies and two reviewers for providing feedback that improved the manuscript. Data analysis and interpretation presented in this study was conducted at the University of Virginia in Charlottesville, Virginia. The University of Virginia was built by enslaved laborers on the unceded lands of the Monacan Nation, who have protected and cultivated these lands for thousands of years. The authors acknowledge and respect their stewardship of the land, past, present, and future.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.863200/full#supplementary-material

References